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addition  of  a  nonlinear  base  isolator.  Finally,  the  method  is  compared  to  MATLAB’s 
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I.  INTRODUCTION 


In  the  wake  of  recent  devastating  earthquakes  in  the  Los  Angeles  and  San 
Francisco  areas,  as  well  as  Mexico  and  Japan,  has  come  heightened  interest  in  earthquake 
protection  systems  for  structures  in  high-risk  areas.  Successful  implementation  of  such 
systems  could  result  in  enormous  savings  in  terms  of  both  dollars  and  lives.  The  goal  is 
not  merely  to  design  a  building  that  can  withstand  the  forces  of  an  earthquake,  but  rather 
to  devise  a  method  to  ensure  that  both  the  building  and  its  contents  survive  intact. 

Simply  strengthening  the  building  would  still  allow  the  vibrational  energy  to  be 
transmitted  to  the  contents,  resulting  in  extensive  internal  damage.  Therefore,  current 
methods  focus  on  systems  that  effectively  isolate  the  entire  building.  These  methods 
consist  of  spring-damper  systems  installed  beneath  the  building  to  absorb  as  much 
earthquake  energy  as  possible. 

Of  critical  importance  is  the  task  of  matching  the  isolator  to  the  structure  in  terms 
of  spring  rate  and  damping  coefficient,  both  of  which  are  generally  nonlinear.  The  design 
of  these  systems  depends  upon  the  mass,  damping,  and  stif&iess  of  the  structure,  all  of 
which  can  be  difficult  to  determine  in  a  complex  structure,  as  well  as  die  frequency  of  the 
groimd  motion.  Performing  actual  vibrational  experiments  on  buildings  is  at  best  cost 
prohibitive,  and,  in  many  cases,  impossible.  In  lieu  of  physical  test  data,  finite  element 
(FE)  techniques  can  be  used  to  construct  and  analyze  mathematical  models  of  structures. 
Following  an  analysis  of  the  structure  model,  the  FE  process  can  then  be  used  to  design 


1 


an  appropriate  isolation  system  and  test  the  response  of  the  modified  structure.  However, 
the  computational  cost  of  conducting  such  an  analysis  can  be  immense  for  a  complex 
structure,  and  this  cost  is  compoimded  when  changing  or  refining  an  isolator  design,  as 
the  entire  process  is  traditionally  repeated  each  time  the  system  is  modified.  Two  current 
methods  of  FE  analysis  include  a  standard  recursive  scheme  and  the  more  recently 
developed  nonlinear  transient  structural  synthesis  method.  Each  of  these  methods  has 
certain  advantages  and  disadvantages  which  will  be  discussed.  Finally,  a  new  method 
will  be  described  which  attempts  to  combine  the  attributes  of  both. 
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n.  EQUATIONS  OF  MOTION  OF  AN  N-STORY  BUILDING 

While  the  methods  discussed  in  this  thesis  can  be  applied  to  a  wide  class  of  finite 
element  problems,  they  will  be  demonstrated  for  analysis  of  building  response  to 
earthquake  excitation.  As  an  example,  the  following  derivation  is  presented  for  a  simple 
four-story  building.  The  procedure  is  easily  extended  to  any  number  of  stories.  A  FE 
model  of  a  large  and  complex  building  would  typically  contain  tens  of  thousands  of  DOF, 
which  is  a  source  of  much  computational  expense. 

In  modeling  buildings,  the  mass  is  commonly  considered  to  be  concentrated  in  the 
floor  of  each  section,  while  the  walls  are  treated  as  massless  colximns  providing  lateral 
stiffiiess  [Ref.  1].  Referring  to  Figure  (1),  the  equations  of  motion  for  each  floor  can  be 
seen  to  be 

OTiii  +(c,  +C2)ii  -C2X2  +(A:i  +A:2)x,  -k2X2  =  /i 
»»2^2  -^2^1  +(^2  +^3)^2  -C3XJ  -k2Xj  +(^2  +k3)X2  -k^X^  =  f 2 
-  C3X2  +  (C3  +  C4  )X3  -  C4X4  -  *3X2  +  (*3  +  *4  )X3  -  *4X4  =  /3 
ffl4X4  -  C4X3  +  C4X4  -  ^4X3  +  A4X4  =  /, 

These  equations  can  be  written  in  matrix  form  as 
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or,  more  compactly, 


(2.1) 


Thus,  for  an  «-story  building,  the  mass  and  stiffiiess  will  be  represented  by  wxn 
matrices.  For  the  sake  of  simplicity  in  demonstrating  the  method,  the  model  has  been 
represented  as  one  degree  of  freedom  per  story.  In  practice,  there  may  be  hundreds  of 
DOF  per  story,  but  the  procedure  remains  the  same.  Ground  motion  is  represented  as 
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y(t).  The  force  transmitted  through  the  isolator  is  generally  nonlinear  and  is  a  function  of 
displacement  and  velocity  across  the  isolator. 
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III.  NONLINEAR  TIME  DOMAIN  STRUCTURAL  SYNTHESIS 

A.  BACKGROUND 

Structiiral  synthesis  refers  to  substructure  coupling  and  structural  modification  in 
a  finite  element  model.  This  discussion  will  provide  an  overview  of  the  structural 
modification  aspect  of  a  recently  developed  structural  synthesis  method,  based  on 
References  [2,3].  Current  work  in  the  area  of  structural  synthesis  centers  on  use  of  a  time 
domain  formulation.  The  goal  of  the  synthesis  method  is  to  reduce  the  computational 
burden  associated  \vith  analyzing  the  effect  of  modifications  to  a  structure.  Without  the 
use  of  synthesis,  the  entire  finite  element  model  must  be  resolved  for  every  variation  in 
the  structure,  involving  potentially  huge  matrix  operations.  With  structural  synthesis,  the 
entire  model  is  solved  only  once,  omitting  any  nonlinearities  (such  as  base  isolators)  in 
the  system.  This  solution  is  relatively  simple  as  it  is  entirely  linear.  The  nonlinearities 
are  accounted  for  as  forces  applied  at  the  c-set  DOF  due  to  relative  displacement  and 
velocity  across  the  isolator.  In  reanalysis,  only  the  cset  DOF  must  be  retained,  using  the 
original  transition  matrices  from  the  baseline  model.  In  this  way,  the  computational 
burden  is  drastically  reduced,  especially  for  successive  solutions  as  in  an  optimization 
routine. 

B.  THEORY 

For  purposes  of  finite  element  analysis,  a  structure’s  physical  coordinates  are 
represented  as  a  vector  {x},  which  is  partitioned  into  {x^}  and  {x^},  with  the  subscript  c 
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referring  to  connection  coordinates,  or  those  coordinates  at  which  the  structure  is  to  be 
modified,  and  the  subscript  i  referring  to  internal  coordinates,  where  no  modification 
takes  place.  In  terms  of  the  convolution  integral,  the  dynamic  response  of  the  system  is 


//,,(/ -r)J[F,(r)J 


(3.1) 


In  structural  modification,  the  interior  coordinates  may  experience  externally 
applied  forces,  while  the  connection  coordinates  may  experience  both  externally  applied 
and  modification  forces.  Thus,  the  force  vector  can  be  represented  as 


(3.2) 


In  Equation  (3.2)  and  throughout  this  discussion,  the  superscript  e  refers  to  externally 
applied  forces,  while  the  superscript  *  denotes  a  quantity  associated  with  the 
modification.  Including  these  synthesized  forces,  the  total  response  becomes 


k(0l  ^ 


which  can  be  rewritten  as 


(3.3) 


k(ol  _|^i(^)L  f 

h 

1 

(3.4) 


The  synthesis  forces  for  linear  structural  modification  can  be  written  generally  as 

Applying  Equation  (3.5)  to  the  second  row  of  Equation  (3.4)  results  in 
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t^c(o}= K(o}- 


(3.6) 

which  is  a  nonstandard  nonhomogeneous  Volterra  integral  equation  of  the  second  kind. 
The  goal  is  to  solve  Equation  (3.6)  for  {xq  (t)},  which  is  the  transient  response  of  the 
modified  system.  In  order  to  obtain  a  numerical  solution,  Equation  (3.6)  can  be 
integrated  by  parts  twice  and  reduced  as  in  Reference  [3],  resulting  in 

j/]+[/(,,(0)lM*fe(o)= 

Equation  (3.7)  can  then  be  solved  numerically  for  {xc*(t)}. 
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IV.  STANDARD  RECURSIVE  METHOD 


Recalling  Equation  (2.1),  we  have  the  governing  differential  equations  for  the 
motion  of  an  «-story  building: 

Multiplying  by  [m"*]  gives 

{x}+([a/-‘Ic])[x}+([m-’1^:])(x}=  (4.1) 

Redefining  {x}  as  the  state  vector  |^| ,  we  have,  in  state-space  notation,  the  differential 


(4.2) 


equation  for  an  n*-order  linear  time-invariant  continuous-time  system: 

{x(0}=[4^(0}+[5if(0} 


M  =  r  r  r  1 


where  x(t)  is  the  n-dimensional  state  vector  and  f(t)  is  the  r-dimetisional  input  vector 
(dropping  brackets  for  clarity).  A  and  B  are  nxn  and  nxr  matrices  of  constant  coefficients, 
respectively.  In  order  to  derive  the  response  of  the  system,  we  multiply  both  sides  of 
Equation  (4.2)  by  the  matrix  K(t),  which  will  be  defined  later: 

mm = mMt) + mBm  (4  3) 

Recognizing  that,  by  applying  the  chain  rule. 


dt 

yfe  see  that 
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^[K{t)x{t)]- k{t)x{t)  =  mAx{t)  +  mBfit) 


dt^ .  (4.4) 

Next,  we  define  K(t)  such  that 

k{t)  =  -Am  (4  5) 

which  has  the  solution 

(4  6) 

We  arbitrarily  choose 

m  =  I  (4.7) 

where  /  represents  the  identity  matrix,  so  that 

m  =  (4.8) 


Recognizing  the  commutability  of  K(t)  and  A,  Equation  4.4  can  be  reduced  to 

^[mxit)]= msm 

at 

Integrating, 


mx{t)  =  ^:(0W0)+  ^K(T)Bf{T)dt 


=  x((!i)+{K{r)Bmdr 

0  (4.10) 

Premultiplying  by  K'^  (t),  we  obtain  the  response 


x{t)  =  ^:-‘(0^(0)  +  K-\t)  \K{r)BmdT 
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=  'I'aWO)  +  J'F(/-r)5/(i-)rfr  (4.1 1) 

0 


where 


(4.12) 


is  known  as  the  transition  matrix.  Thus  the  state  at  a  particular  sampling  time  t  is  given 
by 


xit)  =  e^';c(0)  + 

0 

which  is  a  Volterra  Integro-Differential  Equation  (VIDE)  of  the  second  kind. 
At  the  following  sampling  time,  the  state  is  given  by 

/+A/ 

x{t  +  A/)  = 

0 


(4.13) 


-e 


A/  /+A/ 

e'*‘x(0)+ 


/+A/ 

=  e^^‘{x(,t)}+  (4.14) 

t 

If  the  sampling  period  A/  is  sufficiently  small,  and  the  input  vector  f(t)  is  assumed  to  be 
constant  over  each  time  interval  / :  r  +  A/ ,  the  second  integral  on  the  right  hand  side  can  be 
approximated  as 


i+M  /+A/ 

J  e'^^‘*^-^^BfiT)dT  = 


W(t) 


(4.15) 


Defining  <t  =  ^ + Af  -  t,  the  integral  on  the  right  hand  side  of  Equation  (4. 1 5)  reduces  to 
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So,  from  Equation  (4.14),  we  can  obtain  the  discrete-time  state  vector  sequence 

x(t + AO  =  'FCAOxCO  +  r(  A0/(/)  (4.17) 

where 

'P(A0  =  e""' and  r(A0  =  ^''(e"^ -/> 

Equation  (4.18)  represents  a  recursive  expression  which  will  solve  for  the  state  vector  {x} 
at  each  sampling  time. 

The  above  derivation  is  an  overview  of  that  presented  in  Reference  [4]  and  results 
in  a  recursive  solution  that  has  several  desirable  properties.  IPand  Fzxt  constant 
matrices,  and  so  need  be  computed  only  once  for  any  given  model.  Additionally,  each 
x(t)  and/fO  can  be  discarded  after  ;»:(/+ AO  has  been  computed.  However,  all  DOF  must 
be  retained  for  the  solution.  The  A  matrix  above  can  be  seen  to  be  of  size  2n  x  2n,  where 
n  is  the  number  of  DOF.  Since  formation  of  IPand  /’involve  an  exponential  of  A  and  the 
inverse  of  A,  respectively,  the  computational  requirements  can  be  enormous  for  a  large 
system.  Additionally,  inclusion  of  a  single  nonlinearity  (such  as  a  base  isolator)  renders 


the  entire  model  nonlinear,  and  so  must  be  calculated  as  such. 
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V.  RECURSIVE  SYNTHESIS 


A.  OVERVIEW 

The  goal  of  the  recursive  synthesis  method  developed  in  this  thesis  is  to 
incorporate  the  benefits  of  both  methods  previously  discussed.  Specifically,  it  was 
desired  that  the  method  function  while  retaining  only  the  c-set  DOF,  as  well  as  only 
retaining  the  solution  for  one  previous  time  step.  The  basis  of  the  method  is  a  transition 
matrix  derived  from  the  second-order  differential  equations,  as  opposed  to  the  previously 
used  transition  matrices  based  on  the  homogeneous  solution  to  the  associated  first-order 
differential  equation.  Two  formulations  of  this  method  were  developed,  one  using  a  real 
modal  approach,  the  other  complex. 

Again,  we  begin  with  Equation  (2.1),  the  governing  differential  equation  for  a 
spring  -  mass  system. 

[M]{x}+[c](i}+TO  =  {F} 

Natural  jfrequencies  are  independent  of  damping,  so  we  can  assume  the  following 
solution  to  the  above  equation; 

{x}=y>)C,eJ^  (5.1) 

where  C7  is  an  arbitrary  constant  of  integration  (not  related  to  the  damping  coefficient.) 
Taking  derivatives, 

C,e^  (5.2) 

Substituting  into  Equation  (2.1), 
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-  [M\(f>}6)^C,e^‘^  +  =0 


or 

{K]-a^[M%i,Y:,e“'  =0 


(5.3) 


Now,  we  want  to  solve  the  above  system.  Realizing  that  ^  0  for  finite  t,  and  that  if 
C/  is  zero  we  have  a  trivial  solution,  we  see  that 

The  solution  to  the  above  eigen  system  yields  [»^],  the  diagonal  matrix  of  natural 
frequencies,  and  [i)],  the  modal  transformation  matrix.  [3>]  is  used  to  transform  from 
physical  coordinates  to  modal  coordinates,  as  follows: 


(5.4.a) 

[m]=  [of  [mIcP] 

(5.4.b) 

(5.4.C) 

(5.4.d) 

Transforming  Equation  (2.1)  to  modal  coordinates. 


■\ 

■\ 

\ 

\_ 

(5.5) 


Both  the  real  and  complex  formulations  are  based  on  the  above  modal  differential 
equation. 
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B.  REAL  FORMULATION 


gives 

fe)=[4fe)+{«R  (5.7) 

Although  Equation  (5.7)  is  in  modal  coordinates  and  the  coefficients  are  defined 
differently,  it  is  of  exactly  the  same  form  as  Equation  (4.2)  in  the  derivation  for  the 
standard  recursion.  Therefore,  we  can  proceed  as  before  to  obtain  a  solution  expressed  as 
fe(0}=['f'fej+{r}FW  (5.8) 

where  and  {r}=  ^'^^ ]-[/]){£} 

As  [t,  ]  and  {r,  }  are  constant  for  each  mode,  the  solution  for  the  following  time 
step  is  simply 

{q,(t  +  A/)}  =  rt, )+  {r}fit + At) 
which  leads  us  to  the  recursion: 
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Fit  +  At)  =  Fix(t  +  At)) 


fit  +  At)  =  [<t>YFit  +  At) 


(5.10) 


(5.11) 


(5.12) 


The  key  to  the  synthesis  method  is  in  the  force  term  of  the  recursion.  The  natural 
frequencies  and  transition  matrix  are  determined  from  the  linear  pre-modification  system, 
and  are  retained  when  modifications  are  installed.  The  force  vector  includes  all  forces 
due  to  installation  of  modifications  as  described  in  Chapter  III.  Since  the  recursion  is 
performed  using  the  modal  equations,  which  are  decoupled,  the  analysis  can  be  limited  to 
the  DOF  of  interest  (the  c-set),  providing  substantial  savings  in  computational 
requirements. 

C.  COMPLEX  FORMULATION 


A  complex  formulation  of  the  same  method  has  been  developed  which  has  certain 
computational  advantages  over  the  above  real  formulation. 

Rearranging  Equation  (5.6),  we  can  express  the  modal  equation  of  motion  as 

{qi}+2^0}„.{qi}+(Dl{qi\=\^'}  {^(O},  (5.13) 

We  will  now  find  the  solution  in  the  second  order  form  rather  than  converting  to  a  state 
equation  to  obtain  a  first  order  ODE.  The  total  solution  is  (dropping  brackets) 

^(0  ~  9hom  (0  9 part  (0  (5.14) 

which  for  the  mode  is 
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qi  (0  =  e  [Ai  cos{cdj.  ?)+  £,  sin(0rf.  ^))+  |  hi  (t  -  z)Fi  (r)cfr 


(5.15) 


where; 

B,  =-^fo„ 

A,(0  =  — sin(0^.O 
The  following  are  now  defined: 

#  ±>rf; 


[4(0]= 


0 


1-y 


i+y 


.  1 


{i}=Li  If 


{v}= 


(5.16.a) 

(5.16.b) 


(5.16.C) 

(5.16.d) 


(5.17.a) 


(5.17.b) 


(5.17.C) 


(5.17.d) 

(5.17.e) 


(5.17.f) 
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Using  these  definitions, 


9/(0  =  {ir[4(0l^fe)+  |{ir[i,(?-r)I/’Jv}^(r)/^r 


A,.(o={ir[i,(oio]{v} 

Now,  defining 

fe(o).{-:} 

5 

we  have 

{9/(o}= [i/(oIo]kK  |[i/('-uIoKv/}^(ryr 


(5.18.a) 

(5.18.b) 

(5.19) 
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(5.20) 


which  lead  to  the  total  eqtiation: 

{?(o} = [mlpig, )+  ^[10  -  (5.22) 

Since,  from  Equation  (5.4.d) 

{^(0}=l4>f{F(0}, 
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we  now  have 


We  now  begin  to  develop  the  recursion  with  the  expression  for  the  following  time  step: 

{q{t  +  A/)}  =  [L{t  +  AO][F]{0o  )  +  [l(/  +  A?  -  r)][F][FlcI>]^  {F(r))^/r  (5.24) 

Using  exponential  addition  rules  and  the  commutative  property,  we  see  that 

[!(/,  +/2)]  =  [i(h)I^(^2)]  =  [i(^2)][Ah)]  (5.25) 

so  that 

[q{t  +  A/)}  =  [l(AO][i(0^]{eo  !  +  I"  ^  [i(  {F(r)}rfr 

=  [i(AOli(OM!2o}+  |[i(A0li(f  -  J-)Mf'l4>f  {nr)]dT  + ... 

...  +  1"^  [l(A0ll(/  +  A?  -  T)lFlFl<Df  {F(j)]dT 
=  [I(A/)]J[L(OM0o}+  [{W  -  r)lFlFl<Df  {F(r)Vr|  + ... 

...  +  |"^[Z(A/)Ii(/  +  A?  -  r)lFlFl<I>f  {F(r)}  (5.26) 
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Now,  introducing  the  change  of  variables 

(jst+iu-v,  (5-27) 

we  see  that 

dT  =  -d<7 

and  that  when 

T  =  t,(X=  At 

T  =  t  +  At,a  =  0 

If  At  is  assumed  small,  so  that  {f(0}  is  approximately  constant  over^:/ + A/ ,  we  can  write 
^{t  +  A0}=  [l(AO]|li(/)IP]feo}+  |[^«-  {F(r)]rfr|  +  ... 

...  +  |''[i(cT)]rf<T[/>Ma>f  {f(0}  (5.28) 

which  is 

{^(/+A/)}  =  [i(A/)]{g(o}+[r]{F(/)}  (5.29.a) 

where 

r=  |“[L(c7)]da{pM4.f{F(/)}  (5.29.b) 

So,  the  complete  recursion  is: 


25 


{q{t  +  A/)}  ^  [i(A/)]{9(/)}  +[r]{F(/)} 

(5.30.a) 

{x(r  +  AO}  C=[ci>][i]{^(r+A/)} 

(5.30.b) 

{ F(J  +  A/)}  <=  1  x{t  +  AO} ,  { i(r  +  Ar)} ,  y{t  +  At),  /)| 

(5.30.C) 

t  <=t  +  At 

(5.30.d) 

For  elastic  modes. 

(5.31) 

For  each  term  of  [i,  (cr)]  in  [r],  the  integral  is  easily  evaluated: 

(5.32) 

Substantial  savings  are  realized  by  extracting  from  the  ['P]  matrix  only  those  DOF 
of  interest  (the  cset).  Thus,  the  matrix  becomes  c  x  c  as  opposed  to  2n  x  2n. 

Additionally,  recognizing  that  [i]  is  a  diagonal  matrix,  we  can  further  increase 
computational  efficiency  by  extracting  the  diagonal  terms  and  dot-multiplying  by  the 
terms  of  {?} ,  rather  than  performing  the  full  matrix  multiplication  at  each  time  step. 
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VI.  RESULTS 


Each  formulation  of  this  method  was  compared  to  MATLAB's  ODE45  function,  a 
routine  commonly  used  for  the  solution  of  this  type  of  differential  equation.  Performance 
was  compared  in  the  response  to  a  simulated  blast  loading  produced  by  the  MATLAB 
code  fBlastForcing  .m.  The  forcing  function  is  shown  in  Figure  (2). 


Blast  Force  vs  Time 


Figure  2 

Figure  (3)  shows  the  response  as  computed  by  the  complex  synthesis  method  and 
by  ODE45,  using  40  DOF. 
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Synthesis  vs  ODE45 
(Linear  Spring) 
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Figure  3 

The  plots  of  the  two  solutions  are  indistinguishable  within  the  resolution  of  the 
graph,  demonstrating  the  accuracy  of  the  method.  Plots  using  varying  DOF  as  well  as 
using  the  real  synthesis  method  showed  similar  correlation.  Further  comparisons  were 
made  in  order  to  determine  the  synthesis  method's  efficiency  as  compared  to  ODE45. 
Efficiency  can  be  measured  either  in  terms  of  computational  time  or  number  of  floating 
point  operations  (flops)  required  for  the  solution.  For  both  standards,  a  ratio  was  created 
by  dividing  the  time  (or  number  of  flops)  required  by  ODE45  by  the  time  (or  flops) 
required  by  the  synthesis  method.  Therefore,  a  value  of  unity  on  the  graph  indicates  the 
two  methods  are  equal,  while  any  value  greater  than  one  indicates  a  factor  of 
improvement  by  the  synthesis  method.  In  all  cases,  this  factor  was  plotted  versus  the 
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number  of  degrees  of  freedom  in  the  model.  Figures  (4)  and  (5)  compare  the  real 
formulation  of  the  synthesis  method  using  a  linear  spring  model. 


Linear  Real  Synthesis  Method 
vs  ODE45  (Time) 
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Figure  4 


As  can  be  seen  in  Figure  (4),  ODE45  is  actually  faster  than  the  synthesis  method  for 
models  with  fewer  than  300  degrees  of  freedom.  For  models  of  greater  than  300  DOF, 
the  factor  begins  a  steep  climb,  and  at  500  DOF  the  synthesis  is  slightly  more  than  twice 
as  fast.  As  will  be  discussed  later,  all  time  estimates  are  very  conservative  due  to  the 
adaptive  quadrature  utilized  by  ODE45. 
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Linear  Real  Synthesis  Method 
ODE45  (Operation  Count) 


Figure  5 


Figure  (5)  shows  the  flops  required  for  the  same  comparison  run.  Here,  the  savings  are 
very  noteworthy.  For  the  smallest  model  tested  (four  DOF)  ODE45  required  more  than 
10  times  as  many  flops.  At  500  DOF,  ODE45  required  well  over  600  times  more  flops. 
Again,  even  this  result  is  conservative  due  to  ODE45's  adaptive  quadrature.  Of  great 
importance  is  the  fact  that  in  these  figures,  as  well  as  those  which  follow,  the  factor  of 
improvement  increases  with  the  number  of  DOF.  Thus,  for  a  full-sized  model,  the 
savings  could  potentially  be  several  orders  of  magnitude. 
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The  next  comparison  performed  dealt  with  the  complex  formulation  of  the 


synthesis.  Again,  the  isolator  was  modeled  as  a  linear  spring  with  proportional  damping. 
As  shown  in  Figure  (6),  the  time  savings  are  more  dramatic  than  with  the  real 
formulation.  At  four  DOF,  the  synthesis  is  nearly  ten  times  faster  than  ODE45,  while  at 
400  DOF  it  is  over  115  times  faster. 


Linear  Complex  Synthesis 
Method  vs  ODE45  (Time) 
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Figure  6 

Figure  (7)  shows  that  savings  in  operation  count  are  similar,  with  five  times  fewer 
flops  at  four  DOF,  and  105  times  fewer  at  400  DOF. 
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Linear  Complex  Synthesis  Method 
vs  ODE45 
(Operation  Count) 


Degrees  of  Freedom 


Figure  7 

The  final  trial  again  compared  the  complex  formulation  to  ODE45.  In  this  case, 
the  isolator  was  modeled  as  a  nonlinear  spring  with  proportional  damping  to  more 
accurately  approximate  an  actual  isolator.  The  nonlinearity  was  imposed  using  the 
function  fNonlinearSpring.m,  which  provides  an  interpolated  lookup  table  of  stif&ess 
versus  deflection.  The  force  versus  distance  produced  by  fNonlinearSpring.m  is  shown  in 
Figure  (8).  The  response  obtained  by  the  complex  s5mthesis  method  compared  with 
ODE45  for  40  DOF  is  shown  in  Figure  (9),  demonstrating  its  accuracy.  Again,  the  result 
was  unaffected  by  number  of  DOF  or  by  which  synthesis  formulation  was  used. 
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In  Figure  (10),  the  plot  of  time  factor  versus  DOF,  it  can  be  seen  that  ODE45  is  again 
faster  for  smaller  models,  but  the  synthesis  becomes  faster  at  less  than  100  DOF.  For  400 
DOF,  the  synthesis  is  over  four  times  as  fast. 


Nonlinear  Complex  Synthesis 
Method  vs  ODE45  (Time) 
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Figure  10 


In  terms  of  operational  count,  the  savings  are  again  significant,  with  ODE45  requiring  45 
times  more  flops  at  400  DOF,  as  seen  in  Figure  (11). 
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Figure  11 
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VIL  CONCLUSIONS 


The  second  order  nonlinear  synthesis  method  works,  as  can  be  seen  by  the 
accuracy  of  the  method  compared  to  the  response  generated  by  ODE45. 

Additionally,  the  method  allows  potentially  huge  savings  in  computational 
requirements.  The  complex  formulation  produces  exceptionally  dramatic  time  savings 
due  to  the  use  of  diagonal  matrices.  As  previously  alluded  to,  all  comparisons  in  this 
study  are  very  conservative.  ODE45  uses  adaptive  quadrature  in  solving  the  system  of 
equations.  This  means  that,  if  the  function  is  not  changing  value  significantly  within  each 
time  step,  the  length  of  the  time  step  is  increased,  thereby  reducing  the  total  number  of 
evaluations  required.  The  synthesis  method  as  implemented  in  these  comparisons  does 
not  employ  adaptive  quadrature,  although  the  method  lends  itself  well  to  such  a  scheme. 
With  adaptive  quadrature  installed,  the  synthesis  method  is  expected  to  show  even  more 
dramatic  results  for  all  sizes  of  FE  model. 

Finally,  the  factor  of  improvement  in  all  cases  is  seen  to  increase  widi  increasing 
number  of  DOF.  This  is  compounded  by  the  adaptive  quadrature  issue  discussed  above, 
but  will  still  be  evident  with  adaptive  quadrature  installed.  Due  to  the  reduced  size  of 
matrices  involved  in  the  synthesis,  the  savings  will  be  much  more  pronounced  in  a  FE 
model  of  a  size  more  representative  of  an  actual  building. 
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VIII.  RECOMMENDATIONS  FOR  FUTURE  WORK 


The  real  modal  formulation  of  the  synthesis,  while  extremely  efficient  in  terms  of 
flop  count,  is  less  impressive  in  its  time  savings.  This  is  likely  due  to  the  current 
programming  routine,  which  involves  nested  loops  for  the  solution  of  the  differential 
equation.  If  the  same  method  can  be  programmed  with  fewer  loops,  the  time  savings 
should  improve  commensurately. 

Further  comparisons  should  be  conducted  using  a  more  realistic  representation  of 
the  actual  base  isolators.  In  the  final  trials,  a  nonlinear  spring  was  used  in  conjunction 
with  proportional  damping.  Currently  available  routines  could  be  easily  implemented 
which  model  the  actual  performance  of  base  isolators  very  accurately.  This  would  allow 
definitive  predictions  of  savings  available  through  the  use  of  the  new  synthesis  method  as 
applied  to  earthquake  isolation. 

Finally,  adaptive  quadrature  should  be  included  in  the  synthesis  method.  An 
obvious  next  step,  this  would  be  a  straightforward  modification  to  the  routine,  which 
would  provide  an  even  basis  for  direct  comparisons  with  other  methods. 
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APPENDIX  A:  MATLAB  CODE  FOR  REAL  FORMULATION  OF 


RECURSIVE  SYNTHESIS  WITH  LINEAR  SPRING 


%  realsynth.m  -  Real  formulation  of  recursive  synthesis 
%  method^  using  3-d  matrices.  Uses  linear  spring  for 
%  base  isolator 

%  _ 

clear 

plotme=l ;  compare=l ; 

j  =  sqrt (-1) ; 

global  Amod  B  Yo  k  c  kb  odeforce  tode 

%  TIME  STEPPING: 

%  - 


start_t  =  0.0; 
dt  =0.05; 
end_t  =  20; 

time  =  tstart_t  :dt:end_t]  ; 
n.step  =  length  (time); 


%  Time  points 

%  No.  Time  points 


I  - 

%  Describe  spring-mass  system: 

cset  =  [1] ; 

kel=10000*ones,(l,  500)  ;  %elemental  stiffness 
mel=20000*ones (length (kel) ) ; %elemental  mass 
ndof=length(kel) ; 


k  =  zeros (ndof);  m  =  zeros (ndof) ;  c  =  zeros (ndof); 

%  Populate  [k] ^  [c] , [m] 

k (ndof, ndof)  =  kel (ndof);  m (ndof, ndof)  =  mel (ndof ) ; 
for  i  =  l:ndof-l; 

k(i,i)  =  kel(i)  +  kel(i+l);  k(i,i+l)  =  -kel(i+l); 
k(i+l,i)  =  -kel(i+l) ; 
m(i, i)  =  mel (i) ; 
end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%  Calculate  the  normal  modes  and  natural  frequencies,  and  mass 
normalize 

%  the  eigenvectors.  Sort  the  eigenvalues/vectors  by  ascending 
%  natural  frequency. 

[phi, lam]  =  eig(m\k); 
mtilda  =  phi**m*phi; 
for  i  =  l:ndof 

phi(:,i)  =  phi (:, i) *1/ (sqrt (mtilda (i,  i)  ))  ; 
end 

%  Sort  the  eigenvalues  in  ascending  order. 
ev= (diag (lam) ) ’ ; 

[lam, p] =sort  (ev) ; 
lamstar  =  diag (lam); 
phistar  =  phi; 
for  b  =  l:ndof 
phi ( : , b)  =  phistar ( : , p (b) ) ; 
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end 

Wn  =  sqrt(lam);  %  Radian  natural  frequencies 
freqs  =  Wn/2/pi;  %  Hz 
Zeta  ==  0.0; 

ZetaWn  =  Zeta.^Wn; 

Wd  =  Wn  (1  ~  Zeta.^2)^.5; 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

niomRinodes  =0;  %  number  of  rigid  body  modes 

numEmodes  =  length (lam);  %  number  of  elastic  modes 

phiE  =  phi; 

%  SET  UP  RECURSION: 

%  . - . 

%  ELASTIC  MODES 

%  . . . 


PsiE  =  zeros (2, 2, numEmodes) ;  %  Build  3-dimensional  PsiE  and 
GammaE  matrices 

GammaE  =  zeros (2, 1 , numEmodes) ; 

for  icntEmodes  =  1  :  numEmodes;  %  Will  reference  elastic 

modes  only 

Ae  =  [0  1;- (Wn (icntEmodes) ) ^2  -2*Zeta*Wn(icntEmodes)]; 
PsiE (:,:, icntEmodes )  =  expm(dt*Ae); 

GammaE , icntEmodes ) =Ae\ [PsiE ( : , : , icntEmodes) - 
eye(2)]*[0;l]; 
end 

%  Setup  forcing  function  (base  displacement) : 

Yo  =  1.0;  %  Amplitude 

[y_of_t]  =  fBlastForcing  (Yo,  time  * ,  'blsf,  0); 

%  =  ones (1, nstep) ; 

%  Initialize  force  vector  for  iteration: 

f  =  zeros (l,nstep) ; 
kb  -  15000; 
kc=.l*kb; 

qE  =  zeros (2, 1, numEmodes) ; 

X  =  zeros (2*length(cset),nstep); 


%  RECURSION 

%  . . 

start=flops; 

tic 

for  i-1 : numEmodes 

tempi ( : , i) =GammaE { : , : , i) *phiE (cset , i ) * ; 
end 

for  icnt_tstep  -  1  :  nstep-1; 
for  icntEmodes  =  1  :  numEmodes; 
qE ( : , : , icntEmodes ) =PsiE ( : , : , icntEmodes ) *qE ( ; , : , icntEmodes ) . . . 
ttempl ( : , icntEmodes) *f (icnt^tstep)  ; 
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X  ( : ,  icnt_tstep)  =x  { : ,  icnt_tstep)  +phiE  (cset,  icntEmodes)  *qE  ,  icntEmodes 
); 

end 

f  (icnt_tstep+l)  =”kb*  (x  (1,  icnt_tstep)  ~y_of_t  (icnt_tstep)  )  -.  .  . 
kc*x (2, icnt_tstep) +kel (1) *x (1, icnt_tstep)  ; 
end 

for  icntEmodes  =  1  :  numEmodes; 
qE  ( : ,  : ,  icntEmodes )  =PsiE  { : ,  : ,  icntEmodes )  *qE  { : ,  : ,  icntEmodes )  ,  .  . 

tGammaE  { ; ,  icntEmodes)  *phiE  (cset,  icntEmodes)  ’*f{nstep); 

X  ( : ,  nstep)  “  x  { : ,  nstep)  +phiE  (cset,  icntEmodes)  *qE  { : ,  : ,  icntEmodes)  ; 
end 
toe 

synthflops=flopS“Start 


t  O  *0  O  O  O  1 


ODE45  COMPARISON 


if  compare==l 

start=flops; 

tic 

kchkel=10000*ones (size (kel) ) ;  %elemental  stiffness 
kchkel(l)=kb; 
for  i  =  Irndof-l; 

kchk(i^i)  =  kchkel(i)  +  kchkel(i+l);  kchk(i,i+l)  =  -kchkel (i+l) ; 

kchk{i+l,i)  =  -kchkel (i+1) ; 
end 

kchk (ndof , ndof ) =kchkel (ndof )  ; 
cchkel=zeros (size  (kel) ) ;  %elemental  damping 
cchkel (1) =kc; 
for  i  =  l:ndof-l; 

cchk(i,i)  =  cchkel(i),  +  cchkel(i+l);  cchk(i,i+l)  =  -cchkel  (i+1 )  ; 

cchk{i+l^i)  =  -cchkel (i+1) ; 
end 

cchk (ndof , ndof ) =cchkel (ndof)  ; 

Amod  -  zeros (2*ndof) ; 
odef orce= [ ] ; 
tode=  []; 

7\mod(l:ndof,ndof+l:2*ndof )  =  eye  (ndof); 

Amod (ndof +1 : 2*ndof ^  ndof+1 : 2*ndof )  =  -m\cchk; 

B  =  zeros(2*ndof,ndof); 

B (ndof+1 : 2*ndof, : )  =  inv(m); 

Amod (ndof+1 : 2*ndof , l:ndof)  =  -m\kchk; 

xode  =  zeros (2*ndof, nstep) ; 
xss  =  zeros(2*ndof,nstep); 

[Time,Xode]=ode23  ( ’  vibemdof  ’ ,  time,  xode  ( : ,  1)  )  ; 

odeflops=flops-start 

xode=Xode  * ; 

toe 

end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
if  plotme“l 
figure (1) 

plot (time, X (1, : ) , *  —  * , time, xode (1,  : )  ) 
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legend  ( ’  X_s_y__n_t_h  ’ ,  *  X_0_D_E  ’ ) 
%grid 

title (’ Displacement  vs  Time*) 
xlabel ( *Time  (Seconds) * ) 
ylabel ( *  Displacement ’ ) 
end 
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APPENDIX  B:  MATLAB  CODE  FOR  COMPLEX  FORMULATION  OF 


RECURSIVE  SYNTHESIS  WITH  LINEl^  SPRING 


%  Backdif f 1 .m  -  Complex  formulation  of  recursive  synthesis 
%  method  using  backward  differencing  method  to  obtain  velocity. 
%  Includes  linear  spring  modification 

%  _ 

clear 

j  =  sqrt (-1) ; 

global  Amod  B  Yo  k  c  kb 

plotme==l ;  compare=l ; 

%  TIME  STEPPING: 


start_t  =  0.0; 
dt  =0.05; 
end_t  =  20; 

time  =  [start__t  :dt:end_t]  ’ ;  %  Time  points 

nstep  =  length (time) ;  %  No.  Time  points 

% 

%  Describe  spring-mass  system: 

cset  =  [1]; 

kel=10000*ones (1, 4 ) ;  %elemental  stiffness 
mel=20000*ones (length (kel) ); %elemental  mass 
ndof=length (kel) ; 


k  =  zeros (ndof);  m  =  zeros (ndof);  c  =  zeros (ndof); 

%  Populate  [k] ^ [c] , [m] 

k (ndof, ndof)  =  kel (ndof);  m (ndof, ndof)  =  mel(ndof); 

for  i  =  l:ndof-l; 

k(i,i)  =  kel(i)  +  kel(i+l);  k(i,i+l)  =  -kel(i+l); 
k(i+l,i)  =  -kel(i+l); 
m (i, i)  =  mel  (i) ; 
end 


%  Calculate  the  normal  modes  and  natural  frequencies,  and  mass 
normalize 

%  the  eigenvectors.  Sort  the  eigenvalues/vectors  by  ascending 
%  natural  frequency. 

[phi, lam]  =  eig(m\k); 
mtilda  =  phi**m*phi; 
for  i  =  l:ndof 

phi(:,i)  =  phi (:, i) *1/ (sqrt (mtilda (i, i) )) ; 
end 

%  Sort  the  eigenvalues  in  ascending  order. 
ev= (diag (lam) ) ’ ; 

[lam,p] =sort (ev) ; 
lamstar  =  diag ( lam) ; 
phistar  =  phi; 
for  b  =  l:ndof  ' 
phi ( : , b)  =  phistar ( : , p (b) ) ; 
end 
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Wn  =  sqrt(lam);  %  Radian  natural  frequencies 
freqs  =  Wn/2/pi;  %  Hz 
Zeta  =  0.0; 

ZetaWn  =  Zeta.*Wn; 

Wd  =  Wn  (1  -  Zeta. '"2 ) '^.5; 


SET  UP  RECURSION  FOR  ELASTIC  MODES 


phiE=phi; 

nuniEmodes=ndof ;  %  for  test  structure  -*  no  rigid  body  modes 
P  =  zeros (2*numEmodes, 2*numEmodes) ; 

Ve  =  zeros (2*numEmodes, numEmodes) ; 

Vcol  =  0; 

icntEmodes  =  0; 


for  icntModes  =1:2:  {2*numEmodes-l ) ;  %  Will  reference 

%  elastic  modes  only 

icntEmodes  =  icntEmodes  +1;  %  Start  at  first 

%  elastic  mode 

%  Build  diagonal  vector  of  (non-zero)  eigenvalue  complex 

%  conjugates  (2n  *  1) 

lame  (icntModes)  =  -ZetaWn  (icntEmodes )  +  j  Wd  (icntEmodes)  ; 
lame (icntModes+1)  =  -ZetaWn (icntEmodes )  -  j  *  Wd (icntEmodes) ; 
gammac (icntModes)  =  (exp (lame (icntModes) *dt) -1) /lame (icntModes ) ; 
gammac (icntModes+1) = (exp (lame (icntModes+1 ) *dt) -1) /lame (icntModes+1 ) ; 


%  Build  Tridiagonal  P  matrix  and  Quasi-Block  Diagonal  V  matrix: 

ij  =  [icntModes  icntModes+1];  %  Indices  for  comp  conj  pair 
tempi  =  j *ZetaWn (icntEmodes) /Wd (icntEmodes ) ; 

P(ij,ij)  =  0.5  *  [(1-templ)  (-j /Wd (icntEmodes ));.. . 

(templ+1)  (j /Wd (icntEmodes) ) ]  ; 

Vcol  =  Vcol+1; 

Ve (ij , Vcol)  =  [0;1]; 
one (ij, Vcol)  =  [1;1]; 

end 

Le  =  diag (exp (lamc*dt ) ) ; %  Diagonal  matrix  of  complex  conjugate 
eigenvals 

GammaEx  =  diag (gammac) ;  %  Integral  of  Le  over  dt 
GammaE  =  GammaEx  *  p  *  Ve  *  phiE (cset , : ) ’ ; 


%  Set  up  forcing  function  (base  displacement) : 

Yo  =  1.0;  %  Amplitude 

[y_of__t]  =  fBlastForcing  (Yo,  time  * ,  ’blst’,  0)  ; 

%  Initialize  vectors  for  iteration: 

f  =  zeros (l,nstep) ; 
kb  =  15000; 
kc=.l*kb; 

qE  =  zeros (2*numEmodes, 1) ; %  Non-reduced 
X=zeros (1, nstep) ; 

temp2=phiE (cset , : )  *  one*; 

Le_vector=diag (Le) ; 
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%  RECURSION 

%  - - 

start=flops ; 
tic 

qE  =  Le^vector.*  qE+GainmaE*f  (1)  ; 

X(l)  =  temp2  *  qE; 

Xdot-X(l) /dt; 

f  ( 2 )  =-kb*  {X  ( 1 )  -y_of_t  ( 1 )  )  -kc*Xdot+kel  ( 1 )  *X  ( 1 )  ; 
for  icnt_tstep  =  2  :  nstep-1; 
qE  =  Le_vector.*  qE+GammaE*f(icnt_tstep); 

X  (icnt_tstep)  =  teinp2  *  qE; 

Xdot= (X (icnt_tstep) -X (icnt_tstep-l) ) /dt; 
f (icnt_tstep+l) =-kb* {X (icnt_tstep) -y_of_t (icnt_tstep) )- 
kc*Xdot+kel (1) *X (icnt_tstep) ; 
end 

qE  =  Le  *  qE  +  GammaE  *f(nstep); 

X{nstep)  =  phiE(cset,:)  *  one*  *  qE; 
toe 

synthflops=flops-start 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


ODE45  COMPARISON 


if  compare—l 
start=flops; 
tic 

kchkel=10000*ones {size  (kel) ) ;  %elemental  stiffness 
kchkel(l)=kb; 
for  i  =  l:ndof-l; 

kchk(i,i)  =  kchkel(i)  +  kchkel(i+l);  kchk(i,i+l)  = 

kchk(i+l,i)  =  -kchkel (i+1) ; 
end 

kchk(ndof/ndof ) =kchkel (ndof ) ; 
cchkel=zeros  (size  (kel)  )  ;  %eleinental  damping 
cchkel (l)=kc; 
for  i  =  l:ndof-l; 

cchk(i,i)  =  cchkel (i)  +  cchkel (i+1);  cchk(i,i+l)  = 

cchk(i+l,i)  =  -cchkel (i+1) ; 
end 

cchk (ndof ^ ndof ) =cchkel (ndof) ; 

Amod  =  2eros(2*ndof); 

Amod(l:ndof,ndof+l:2*ndof )  =  eye(ndof); 

Amod {ndof+l:2*ndof, ndof +l:2*ndof)  =  -m\cchk; 

B  =  2eros(2*ndof,ndof); 

B (ndof+l : 2^ndof , : )  =  inv(m); 

Amod (ndof +1 : 2*ndof, 1 :ndof)  =  -m\kchk; 
xode  =  zeros (2*ndof, nstep) ; 

[Time,Xode] =ode23 { *  fvibemdof * , time,  xode ( : ,  1)  )  ; 

odeflops=flops-start 

xode=Xode  * ; 

toe 

end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
if  plotme==l 


-kchkel (i+1 ) 


-cchkel (i+l) 
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figure (1) 

plot  (tiine,X,  ’  —  \  time,  xode  (1,  : )  ) 
legend  ( *  X_s_y_n_t_h  * ,  *  X_0_D__E  * ) 
title (’ Displacement  vs  Time*) 
xlabel { *  Time  (Seconds ) * ) 
ylabel ( ’ Displacement  * ) 
end 
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APPENDIX  C:  MATLAB  CODE  FOR  COMPLEX  FORMULATION  OF 


RECURSIVE  SYNTHESIS  WITH  NONLINEAR  SPRING 


%  Backdiff2.in  -  Complex  formulation  of  recursive  Synthesis 
%  method  using  backward  differencing  to  obtain  velocity 
%  fNonlinearspring  used  for  modification 

Q, 

"O 

%  _ 

clear 

j  =  sqrt (-1) ; 

global  Amod  B  Yo  k  c  kb 

plotme=l;  compare=l; 

%  TIME  STEPPING: 

%  - 

start__t  ==  0.0; 
dt  =  0 . 0  5 ; 
end__t  =  20; 

time  =  [start_t  :dt  :end_t]  ’ ;  %  Time  points 

nstep  =  length(time) ;  %  No.  Time  points 


%  Describe  spring-mass  system: 

cset  =,[!]; 

kel=10000*ones (1, 4) ;  %elemental  stiffness 
mel=20000*ones  (length (kel) ) ; %elemental  mass 
ndof =length ( kel ) ; 

%  - 

k  =  zeros (ndof);  m  =  zeros (ndof);  c  =  zeros (ndof); 

%  Populate  [k] , [c] ,  [m] 

k(ndof,ndof)  =  kel(ndof);  m(ndof^ndof)  =mel(ndof); 

for  i  =  l:ndof-l; 

k(i,i)  =  kel(i)  +  kel(i+l);  k(i,i+l)  =  -kel(i+l); 
k(i+l,i)  =  -kel(i+l); 
m(i, i)  =  mel (i) ; 
end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%  Calculate  the  normal  modes  and  natural  frequencies,  and  mass 
normalize 

%  the  eigenvectors.  Sort  the  eigenvalues/vectors  by  ascending 
%  natural  frequency. 

[phi, lam]  =  eig(m\k); 
mtilda  =  phi’*m*phi; 
for  i  =  l:ndof 

phi(:,i)  =  phi (:,i)*l/ (sqrt (mtilda (i, i)  ))  ; 
end 

%  Sort  the  eigenvalues  in  ascending  order. 
ev= (diag (lam) ) ' ; 

[lam, p] =sort (ev) ; 
lamstar  =  diag ( lam) ; 
phistar  =  phi; 
for  b  =  l:ndof 
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phi(:,b)  =  phistar(:,p(b) ) ; 
end 

Wn  =  sqrt(lam);  %  Radian  natural  frequencies 
freqs  =  Wn/2/pi;  %  Hz 
Zeta  =  0.0; 

ZetaWn  =  Zeta.*Wn; 

Wd  =  Wn  .*  (1  -  Zeta. "^2)^.5; 


SET  UP  RECURSION; 


% 

%  ELASTIC  MODES 


phiE=phi; 
numEraodes=ndof ;  % 
P 

Ve 

Vcol 

icntEmodes 


for  test  structure  -  no  rigid  body  modes 
=  zeros {2*numEmodes, 2*numEmodes) ; 

-  zeros {2*numEmodes, numEmodes ) ; 

=  0; 

=  0; 


for  icntModes  =  1:2:  (2*n;amEmodes“l )  ;  %  Will  reference 

%  elastic  modes  only 

icntEmodes  =  icntEmodes  +1;  %  Start  at  first 

%  elastic  mode 

%  Build  diagonal  vector  of  (non-zero)  eigenvalue  complex 

%  conjugates  {2n  *  1) 

lame (icntModes)  =  -ZetaWn (icntEmodes)  +  j  *  Wd (icntEmodes) ; 

lame (icntModes+1)  =  -ZetaWn (icntEmodes)  -  j  *  Wd (icntEmodes) ; 
gammac (icntModes)  =  (exp (lame (icntModes) *dt) -1) /lame (icntModes) ; 
gammac (icntModes+1 ) = (exp (lame (icntModes+1) *dt ) - 
1) /lame (icntModes+1) ; 


%  Build  Tridiagonal  P  matrix  and  Quasi-Block  Diagonal  V 

matrix: 

ij  =  [icntModes  icntModes+1] ;  %  Indices  for  comp  conj  pair 
tempi  =  j * ZetaWn (icntEmodes) /Wd (icntEmodes) ; 

P(ij,ij)  =  0.5  *  [ (1-templ)  (- j /Wd (icntEmodes) );.. . 

(templ+1)  (j/Wd( icntEmodes) ) ] ; 

Vcol  =  Vcol+1; 

Ve(ij,Vcol)  =  [0;1]; 
one (ij, Vcol)  =  [1;1]; 

end 

Le  =  diag (exp (lamc*dt) ) ;%  Diagonal  matrix  of  complex  conjugate 
eigenvals 

GammaEx  =  diag (gammac) ;  %  Integral  of  Le  over  dt 
GammaE  =  GammaEx  *  p  *  Ve  *  phiE (cset, : )  * ; 


%  Set  up  forcing  function  (base  displacement) : 

Yo  =  1.0;  %  Amplitude 

[y_of_t]  =  fBlastForcing (Yo, time  * ,  *blst’,  0); 

%  Initialize  vectors  for  iteration: 

f  =  zeros (l,nstep) ; 
kb  =  15000; 


50 


kc=0*kb; 

qE  =  zeros  (2*nuinEmodes,  1)  ;  %  Non-reduced 
X=2eros (l,nstep) ; 
temp2=phiE (cset, : )  *  one'  ; 
Le_vector=diag (Le) ; 

%  RECURSION 


start=flops; 

tic 

qE  =  Le_vector.*  qE+GammaE^f ( 1) ; 

X{1)  =  temp2  *  qE; 

Xdot=X(l) /dt; 

f (2) =-fNonlinearSpring (X{1) -y_of_t (1) , 0) ~kc*Xdot+kel (1) *X (1) ; 
for  icnt_tstep  =  2  :  nstep-1; 
qE  =  Le_vector.*  qE+GainmaE*f  (icnt_tstep)  ; 

X (icnt_tstep)  =  temp2  *  qE; 

Xdot= {X (icnt^tstep) -X (icnt_tstep~l) ) /dt ; 
f  (icnt_tstep+l)  =-fNonlinearSpring  (X  (icnt^tstep) 
y__of_t  (icnt_tstep)  ,  0)  -  kc*Xdot+kel  (1)  *X  (icnt_tstep)  ; 
end 

qE  =  Le  *  qE  +  GammaE  *f(nstep); 

X(nstep)  =  phiE(cset,:)  *  one’  *  qE; 
toe 

synthflops=flops-start 


ODE45  COMPARISON 


j 


if  compare==l 
start-flops; 
tic 

kchkel=10000*ones  (size  (kel)  ) ;  %elertiental  stiffness 
kchkel (1) =0; 
for  i  =  l:ndof-“l; 

kchk{i,i)  =  kchkel (i)  +  kchkel (i+1);  kchk(i,i+l)  =  -kchkel (i+1) 

kchk(i+l,i)  =  -kchkel (i+1) ; 
end 

kchk (ndof , ndof ) =kchkel (ndof ) ; 
cchkel=zeros (size (kel) ) ;  %elemental  damping 
cchkel (1) =kc; 
for  i  =  l:ndof-l; 

cchk(i,i)  =  cchkel (i)  +  cchkel (i+1);  cchk(i,i+l)  =  -cchkel (i+1) 

cchk(i+l,i)  =  -cchkel (i+l) ; 
end 

cchk (ndof , ndof) =cchkel (ndof)  ; 

Amod  =  zeros(2*ndof); 

Amod(l:ndof,ndof+l:2*ndof )  =  eye (ndof); 

Amod (ndof +1 : 2*ndof , ndof +1 : 2*ndof )  =  -mXcchk; 

Amod (ndof +1 : 2*ndof, 1: ndof)  =  -m\kchk; 

B  =  zeros (2*ndof, ndof ) ; 

B (ndof+1 : 2*ndof , : )  =  inv(m); 
xode  =  zeros (2*ndof, nstep) ; 

[Time,Xode] =ode45 ( ’ fvibemdof2 ’ ^  time, xode (:,1)); 

odeflops=flops-start 

xode=Xode ’ ; 
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toe 

end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
if  plotme==l 
figure (1) 

plot (time,X, '  —  ' , time,  xode (1,  ; ) ) 
legend { ' X_s_y_n_t_h ' , ' X_0_D_E ' ) 
title (' Displacement  vs  Time') 
xlabel('Time  (Seconds)') 
ylabel ( ' Displacement ' ) 
end 
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APPENDIX  D:  FUNCTIONS  CALLED  BY  ABOVE  CODES 


function  [f_of_t, fdot]  =  fBlastForcing (Fo, time,  type,  plotit); 

□ 


□ 

%  Usage:  [f_of_t, fdot]  =  fBlastForcing (Fo, time,  type,  plotit); 

□ 

% 

□ 

%  Choices:  sine  blst  step 

□ 

% 

□ 

%  type  =  'step'  STRING  Variable 

□ 

%  - 

□ 

% 

□ 

% 

□ 

%  If  use  'sine*,  fdot  also  returned. 

% 

%  This  function  returns  a  forcing  function  which  is 
%  a  "blast"  function. 

% 

%  F(t)  =  Fo  *  (  exp (-at)  -  exp(-bt)  ) 

% 

%  where  a  and  b  are  constants  which  shape  the  blast, 

%  and  Fo  is  the  amplitude  of  the  blast. 

% 

%  The  variable  "plotit"  is  a  switch  which  if  set  to  an 

%  integer  greater  than  0  will  cause  f(t)  to  be  plotted, 

%  in  the  figure  window  with  that  number,  i.e.  figure (plotit ) . 

%  if  set  to  anything  else,  will  not  plot. 

% 

% 


%  Choices:  sine  blst  step 

%  type  =  'step'; 

%  - 

if  type  —  'blst*; 

disp(*  Blast  forcing  used...*) 
a  =  2.0; 
b  -  5.0; 

f_of_t  =  Fo  *  (  exp(-a*time)  -  exp(-b*time)  ); 

fdot  =  Fo  *  (-a*  exp(-a*time)  +  b*exp (-b*time)  ); 

elseif  type  —  'step'; 
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disp('  Step  forcing 


used. . . * ) 


f_of_t  =  Fo  *  ones  (size  (tiine)  )  ; 
fdot  =  []; 

elseif  type  ==  ’sine*; 

disp(’  Sine  forcing  used...’) 

W  =  1;  %  Hertz 

disp (sprintf ( ’  Freq  (Hz):  %5.1f’,W)) 
f_of_t  =  Fo  *  sin (2*pi*W*time) ; 

fdot  -  Fo  *  (2*pi*W) *cos (2*pi*W*time) 


end; 

if  plotit  >  0; 

figure (plotit) 

if  type  ==  ’sine’; 

plot (time, f_of_t, time, fdot) ;grid 

elseif  type  ===  ’blst’; 

plot (time, f_of_t, time, fdot) ;grid 

else 

plot  (time,  f_of_t)  ;grid 
end 

end 

%  End  function. 
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function [ f ofx] =fNonlinearSpring (xin, klinf it,  plotit ) 


skip='y' ; 
if  skip~='y' 

xvsf= [0 .  C 

)  0 

0.1 

100.0; 

0.2 

200.0; 

0.3 

300.0; 

0.4 

400.0; 

0.5 

500.0; 

0.6 

580.0; 

0.7 

640.0; 

0.8 

700.0; 

0.9 

740.0; 

1.0 

760.0; 

1.1 

780.0; 

1.2 

790.0; 

1.3 

800.0; 

1.4 

805.0; 

1.5 

807.5; 

1.6 

808.75; 

1.8 

809.0; 

10.0 

850.0; 

100.0 

1000.0; 

190.0 

1150.0; 

280.0 

1300,0; 

370.0 

1450.0; 

460.0 

1600.0; 

550.0 

1750.0; 

640.0 

1900.0] 

end 

xvsf= [0 . 0 

0.0; 

0.1 

200.0; 

0.2 

400.0; 

0.3 

600.0; 

0.4 

780.0; 

0.5 

960.0; 

0.6 

1140.0; 

0.7 

1320.0; 

0.8 

1500.0; 

0.9 

1600.0; 

1.0 

1650.0; 

1.1 

1700.0; 

1.2 

1740.0; 

1.3 

1750.0; 

1.4 

1760.0; 

1.5 

1765.0; 

1.6 

1770.0; 

1.8 

1770.0; 

10.0 

1770.0; 

100.0 

1770.0; 

1000.0 

1780.0; 

1900.0 

1790.0; 

2800.0 

1800.0; 

3700.0 

1810.0; 

4600.0 

1820.0; 
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5500.0  1830.0; 

6400.0  1840.0; 

7300.0  1850.0; 

8200.0  1860.0; 

9100.0  1870.0; 

10000.0  1880.0]; 

num_pts=length (xvsf ) ; 

[x, f ] =fXYreflect (xvsf ( : , 1) , xvsf ( : ,  2)  )  ; 
if  klinfit>0; 

f linf it=klinf it*xvsf ( : ,  1)  ; 

[x, flinfit] =fXYreflect (xvsf ( : , 1) , flinfit )  ; 
fdif=f-f  Unfit ; 
else 

fdif=zeros (size (f ) ) ’ ; 
flinfit^zeros (size (f ) ) * ; 

end 

f ofx=interpl (x, f , xin) ; 
if  nargin==3; 

plot (x, f,  ’ r *  r  X, flinfit, * g* , x, fdif , ’y’ ) ;grid; figure (gcf ) 
title ( 'Nonlinear  Spring (red) -Linear  Fit  Spring (grn) - 
Difference (yel) ' ) 

xlabel { ' Deflection  (in) ’ ) 
ylabelC  Force  (Ibf)  ») 

end 
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function  [xreflect, yreflect]  =  fXYreflect (x, y) ; 

%  fXYreflect .m 

%  -  % 

%  Usage:  [xreflect , yreflect]  =  fXYreflect (x, y) ;  % 

%  This  function  creates  a  new  set  of  x,y  data  from  a  given  set  of  x,y 
data. 

%  The  vectors  x  and  y  are  inputted,  and  these  data  are  reflected  about 
%  the  y  axis.  The  new  data  xreflect, yxreflect  contains  both  the 
original 

%  and  the  reflected  data,  and  hence  the  new  vectors  are  of  length 
%  2  *  length (x)  -  1 

%  The  function  requires  that  the  x  vector  start  at  zero 
%  length (x)  ==  length (y)  % 

% 


if  length (x)  ==  length (y)  &  x(l)  ==  0; 
num_pts  =  length(x); 
xreflect (l:num_pts)  =  -flipud(x); 
xreflect  (nuni_pts+l:2*num_ptS“l)  =  x{2:num_pts); 
yreflect (l:num_pts)  =  -flipudCy); 
yreflect (num_pts+l : 2*num_ptS“l)  =  y (2 :num_pts) ; 

end 
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%fvibemdof .m 

%  Used  for  ODE  comparison  -  Linear  spring  modification 
function  xdot=fvibemdof (t , x) 
global  Amod  B  Yo  k  c  kb 
ndof=length (Amod) /2; 

Force=Eeros (ndof , 1) ; 

b  =  5.0; 

a=2.0; 

F_of__t=Yo*  (exp  (~a*t )  -exp  (-b*t )  )  ; 
fdot  =  Yo* (-a*exp(-a*t)+b*exp(-b*t) ) ; 

Force  (1)  =  kb*F__of_t; 
xdot=zeros (size (x) ) ; 
xdot=Amod*x+B* Force; 
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%  fvibeindof2  .m 

%  Used  for  ODE  comparison  with  fNonlinearSpring 
function  xdot=fvibemdof2 (t ,  x) 
global  Amod  B  Yo  k  c  kb 
ndof=length (Amod) /2; 

Force=zeros (ndof , 1) ; 
b  =  5.0; 
a=2 . 0 ; 

F___of__t=Yo*  (exp  (-a*t)  -exp  (~b*t)  )  ; 
fdot  ==  Yo*  {-a*exp (-a*t) +b*exp  (-b*t )  )  ; 

Force(l)  =  -fNonlinearSpring (x {1 ) -F_of_t,  0) ; 

xdot=zeros (size (x) ) ; 

xdot=Amod*x+B*Force; 
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